﻿#pragma once
#include <vtkm/cont/ArrayHandle.h>         //数组类的头文件 头文件是路径
#include <vtkm/worklet/WorkletMapField.h>  //算子封装头文件
#include <vtkm/cont/Algorithm.h>
#include <vtkm/Matrix.h>
#include <iostream>

//有了“clang-format”文件，Ctrl+j+d 自动排版
//using vtkm::FloatDefault; //using用法！！！！！
using Real = vtkm::FloatDefault;
using vtkm::Id;


struct Mxpy :vtkm::worklet::WorkletMapField //一一映射的工作组
{
	using ControlSignature = void(FieldIn x, FieldInOut y);
	using ExecutionSignature = void(_1, _2);                   // 这里输出的值为数组内，某一个下标号里面的值
    
	Real _alpha;
	Mxpy(const Real& alpha) :_alpha(alpha) {}

	template<typename FieldType>
	void operator()(const FieldType& x, FieldType& y) const
	{
		y += x * _alpha;
	}
};

template<typename ArrayType> 
void mxpy(const ArrayType& x, ArrayType& y, const Real& alpha)
{
	vtkm::cont::Invoker{}(Mxpy{ alpha }, x, y); //这里的Invoker调用的是前面的operator！！！！！？
}

template<typename ArrayType>
void fill(ArrayType& x, typename ArrayType::ValueType value)
{
	vtkm::cont::Algorithm::Fill(x, value);         // 这里输出的值为数组内，某一个下标号里面的值（如 x=8，value=1时，表示x数组中x[8]=1）！！！
}

/*
struct Axpby :vtkm::worklet::WorkletMapField
{
	using ControlSignature = void(FieldIn x, FieldIn y, FieldOut z);                            
	using ExecutionSignature = void(_1, _2, _3);                                                      //这里有一个括号赋值法，“：_alpha(alpha), _beta(beta)"
	FloatDefault _alpha, _beta;                                                                       // _alpha:下标在这里是表示:这里单纯是变量定义，并不是固定
	Axpby(const FloatDefault& alpha, const FloatDefault& beta) :_alpha(alpha), _beta(beta) {}  //这个构造函数实在给标量赋值       // 这里axpby是调用的 还是重新定义的(这里的{})：这里是结构体下面的构造函数，要与后面的函数命名区别开来，Axpby和axpby
	                                                                                             //正常的话"_alpha(alpha), _beta(beta)"应写在{}里面，这个简单写在外面属于优化
	template<typename FieldType>
	void operator()(const FieldType& x, const FieldType& y, FieldType& z) const                          
	{
		z = _alpha * x + _beta * y;
	}
};
template<typename ArrayType>
void axpby(const ArrayType& x, const ArrayType& y, ArrayType& z, const FloatDefault& alpha, const FloatDefault beta)
{
	vtkm::cont::Invoker{}(Axpby{alpha,beta}, x, y, z); //vtkm::cont::Invoker{}帮助文件简写，Invoker调用时第一个参数必须是worklet里的，后面的多个参数变量
	                                                   //按需提取（通常是：ControlSignature）
}


#include<vtkm/cont/ArrayHandle.h>
#include<vtkm/worklet/WorkletMapField.h>
#include<vtkm/cont/Algorithm.h>

using vtkm::FloatDefault;
struct AXPBYPCZ:vtkm::worklet::WorkletMapField
{
	using ControlSignature=void(FieldIn x,FieldIn y,FieldIn z,FieldOut W);
	using ExecutionSignature=void(_1,_2,_3,_4);
    FloatDefault _alpha, _beta, _gamma;
    AXPBYPCZ(const FloatDefault& alpha,const FloatDefault& beta,const FloatDefault& gamma):_alpha(alpha), _beta(beta), _gamma(gamma){}

	template<typename FieldType>
    void operator()(const FieldType& x,const FieldType& y,const FieldType& z, FieldType w) const
    {
          w = alpha * x + beta * y + gamma * z;
    }
};
template<typename ArrayType>
void axpbypcz(const ArrayType& x,const ArrayType& y,const ArrayType& z,ArrayType w,const FloatDefault& alpha,\
	const FloatDefault& beta,const FloatDefault& gamma)
{
  vtkm::cont::Invoker{}(AXPBYPC{ alpha, beta, gamma }, x, y, z, w);
}
*/

/*
稀疏矩阵
*/

//#define N 6
//#define M 6
/*
using vtkm::FloatDefault;
struct SparseMatrix : vtkm::worklet::WorkletMapField
{
  using ControlSignature = void(FieldIn& A, FieldOut& B, FieldOut& C, FieldOut& D);
  using ExecutionSignature = void(_1, _2, _3, _4); //用于定义输入矩阵和输出的三个向量组

  FloatDefault _r, _c, _v; //行列参数类型定义
  SparseMatrix(FloatDefault& r, FloatDefault& c, FloatDefault& v): _r(r) , _c(c) , _v(v){} //用来定义 行 列 元素值
  template<typename ArrayType>
  void operator()(const ArrayType& A, ArrayType& B,ArrayType& C, ArrayType& D,FloatDefault& r,FloatDefault& c,FloatDefault& v) const
  {
    //  ①引入稀疏矩阵

	  vtkm::Matrix<FloatDefault, 4, 4> A;
          A[0][0] = 0, A[0][1] = 0, A[0][2] = 0, A[0][3] = 1;
          A[1][0] = 0, A[1][1] = 2, A[1][2] = 0, A[1][3] = 0;
          A[2][0] = 0, A[2][1] = 0, A[2][2] = 1, A[2][3] = 0;
          A[3][0] = 4, A[3][1] = 0, A[3][2] = 0, A[3][3] = 5;
   /*   
   创建了下列的稀疏矩阵A：
          { 0, 0, 0, 1,},
          { 0, 2, 0, 0 },
          { 0, 0, 1, 0 },
          { 4, 0, 0, 5 }, 
   */
    
   // ② 遍历稀疏矩阵 找出非0元素

